function [ExoTimeStates,RegionsForExoStates,pproc] = GenExoTimeStates(DS,opt,PriceChange)
% Create the time varying exogenous state variables indices
% pproc is the key 

% Load price data:

load([opt.DataFolder opt.PricesData])

% Set ARIMA model for sugar price:
MdlPS = arima(1,0,0);

% Get sugar price used:
sug = double(prices(:,opt.SugarPrice));

% Estimate AR(1) for sugar:

%disp(' Estimate process for sugar prices: ')
MdlPS = estimate(MdlPS,sug,'Display','off');

if PriceChange.on == 0
    % Price process:
    pproc.sugar.Cte    = MdlPS.Constant; % This is the the constant.
    pproc.sugar.AR     = MdlPS.AR{1};
    pproc.sugar.Sigma  = sqrt(MdlPS.Variance);
elseif PriceChange.on == 1 % This will change the process for prices, this will be
                           % useful for counterfactual
    pproc.sugar.Cte    = MdlPS.Constant*PriceChange.factor;
    pproc.sugar.AR     = MdlPS.AR{1};
    pproc.sugar.Sigma  = sqrt(MdlPS.Variance);
end    

% Build outside option short series (this will be used to match actual decision):
% Select prices for years used in the estimation:
prices_used = [];
for y = [2003:2013] % Those are all available years of canasat data. opt.Years will be used later for selection in estimation
    prices_used = [prices_used; double(prices(prices.year == y, {opt.CornPrice opt.SoyPrice opt.SugarPrice}))];
end

% Get results from outside option model estimation:
if opt.WhichOutside == 1
    [OutDS, OutOptionGrid] = OutsideOptionModel(DS,opt);
elseif opt.WhichOutside == 2
    [OutDS, OutOptionGrid] = OutsideOptionModelSimple(DS,opt);
elseif opt.WhichOutside == 3
    [OutDS, OutOptionGrid] = OutsideOptionModelOld(DS,opt);
end

OutIndices = unique(OutDS.Index)';

ExoTimeStates        = nan(size(DS,1),size(prices_used,1));
RegionsForExoStates  = nan(size(DS,1),1);

for r = OutIndices
    
    if r ~= 0
    % Build outside option long series (this will be used in the estimation of the outside process):
    % out_return has the shape: r(p,h) = r(p_ave,h) + dr(p_ave,h)/dp * (p - p_ave)
    out_return(:,r) = OutOptionGrid.Base(r) + OutOptionGrid.Partialcorn(r)*(double(prices(:,opt.CornPrice)) - OutOptionGrid.PriceAverage(1)) ...
        + OutOptionGrid.Partialsoy(r)*(double(prices(:,opt.SoyPrice)) - OutOptionGrid.PriceAverage(2));
    
    % Work with return in (R$ 000):
    out_return(:,r) = out_return(:,r)/10^3;
    
    % Create OuReturn variable to be matched to data:
    OutReturn(:,r) = OutOptionGrid.Base(r) + OutOptionGrid.Partialcorn(r)*(prices_used(:,1) - OutOptionGrid.PriceAverage(1)) ...
        + OutOptionGrid.Partialsoy(r)*(prices_used(:,2) - OutOptionGrid.PriceAverage(2));
    OutReturn(:,r) = OutReturn(:,r)/1000;
    
%    disp([' Estimate process for class ' num2str(r) ' return: '])
    % Estimate AR(1) for outside option return:
    MdlOR = arima(1,0,0);
    MdlOR = estimate(MdlOR,out_return(:,r),'Display','full');
    
    % Price process:
    pproc.outside{r}.OutIndex = r;
    pproc.outside{r}.Cte      = MdlOR.Constant; % This is the the constant.
    pproc.outside{r}.AR       = MdlOR.AR{1};
    pproc.outside{r}.Sigma    = sqrt(MdlOR.Variance);
%    disp(' ')
%    disp(' ')
    
    % Discretize the price process by region:
    
    % The process for sugar price and return of outside option is given by
    % the following VAR(1):
    
    % y_t = F0 + F * y_{t-1} + \varepsilon
    % \varepsilon \sim Sigma
    
    Sigma = [pproc.sugar.Sigma^2,  0; 0, pproc.outside{r}.Sigma^2];
    F     = [pproc.sugar.AR, 0; 0, pproc.outside{r}.AR];
    F0    = [pproc.sugar.Cte; pproc.outside{r}.Cte];
    
    % Discretize process:
    [pproc.P{r}, pproc.y{r}, pproc.x{r}, MVndx] = VarTauchen(opt.pproc.Ns,F0,F,opt.pproc.bandwidth,Sigma);
    
    % Now we need to classify prices in the data:
    SugOut = [prices_used(:,3), OutReturn(:,r)];
    
    Ny = size(pproc.y{r},1);
    Nd = size(SugOut,1);
    
    K1 = kron(SugOut(:,1),ones(1,Ny)) - kron(pproc.y{r}(:,1)',ones(Nd,1));
    K1 = abs(K1);
    
    K2 = kron(SugOut(:,2),ones(1,Ny)) - kron(pproc.y{r}(:,2)',ones(Nd,1));
    K2 = abs(K2);
    
    
    [~,exo_idx(:,r)] = min(K1+K2,[],2); 
    clear K1 K2
    
    NumInOut = sum(OutDS.Index == r);
    ExoTimeStates(OutDS.Index == r,:) = repmat(exo_idx(:,r)', [NumInOut,1]);
    RegionsForExoStates(OutDS.Index == r) = repmat(r, [NumInOut,1]);
    
    end
    
    
end


%% plot price of sugar and return of outside option (Article Figure 4):

figure(1)
plot(prices.year(36:end),prices.sb1brl_lag(36:end),'LineWidth',2,'Color','black')
grid on
ylim([0 .5])
%title('Sugar price - Sugar #11')
xlabel('Years')
ylabel('R$/lb.')
set(gca,'fontname','times', 'FontSize', 11);
saveas(gcf,'Figures/SugarPrice.eps', 'epsc')
figure(2)
plot(prices.year(37:end),1000*out_return(37:end,:),'LineWidth',2,'Color','black')
%title('Other uses return - Different estimation bins')
grid on
xlabel('Years')
ylabel('R$/ha.')
set(gca,'fontname','times', 'FontSize', 11);
saveas(gcf,'Figures/OutReturn.eps', 'epsc')

end




